snakes%>%
  group_by(Species)%>%
  count(name = "N")

snakes%>%
  filter(Gear.Type != "")%>%
  select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
  distinct()%>%
  # Calculating tow duration
  mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
                                        is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
  group_by(Gear.Type)%>%
  summarise(n = n(), # counting number of trips sampled
            haul.hours = sum(Tow.duration.hours, na.rm = T)) # total fishing effort

# calculating the mean SVL across years

age.yr <- snakes%>%
  group_by(Species, Year)%>%
  summarise(N = n(),
            mean = mean(Snout.to.Vent..cm., na.rm = T))

# SVL cut off for age classes

maturtity <- snakes%>%
  group_by(Species)%>%
  count()%>%
  mutate(juv = 35,
         adult = ifelse(Species == "Hydrophis curtus", 54, 65))

snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  #calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

# percentage of gravid females in sample

HS%>%
  count(Gravid)%>%
  mutate(N = sum(n))%>%
  mutate(prop.gravid = n/N)%>%
  filter(Gravid == "Yes")%>%
  select(-c(Gravid))

# checking the number of gravid females per year

HS%>%
  group_by(Year)%>%
  filter(Gravid == "Yes")%>%
  count(Gravid)%>%
  spread(Gravid, n)%>%
  rename(`n(Gravid Females)` = Yes)

# calculating the proportion of gravid females per month

gravid <- HS%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
  group_by(Month)%>%
  summarise(N = n(),
            gravid = sum(Gravid == "Yes"),
            prop.gravid = gravid/N)

gravid

# plotting prop gravid per month

gravid%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month, fill = list(prop.gravid = 0))%>%
  ggplot(aes(Month, prop.gravid))+
  geom_point(size = 3)+
  geom_path(aes(group = 1), size = 1, linetype = "dashed")+
  geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Proportion of gravid females")

embryos%>%
  summarise(n.mothers = length(unique(Field.Code)),
            n.embryos = length(unique(Embryo.Code)))

embryos%>%
  select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
  skimr::skim()%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

embryos%>%
  group_by(Field.Code)%>%
  filter(Sex != "")%>%
  summarise(prop.female = sum(Sex == "Female")/n())%>%
  skimr::skim(prop.female)%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

embryos%>%
  count(Embryo)%>%
  mutate(N = sum(n), rate = n*100/N)%>%
  select(-N)

re_clutch <- re%>%
  select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
  group_by(Field.Code)%>%
  summarise(clutch.size = n(),
            Clutch.wt = sum(Egg.Weigth..g..),
            Mother.wt = last(Mother.Wt),
            Female.SVL = last(Female.SVL))%>%
  mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))

re_clutch%>%
  ggplot(aes(Female.SVL, clutch.size))+
  geom_point(size =3)+
  geom_smooth(method = 'glm', method.args = list(family = 'poisson'))+
  labs(x = "Female SVL (cm)", y = "Clutch Size")

mcfadden <- function(x){
  
  r2 <- 1 - (x$deviance/x$null.deviance)
  
}

re_clutch%>%
  select(Female.SVL, clutch.size)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(clutch.size ~ Female.SVL, data = ., family = "poisson")),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance),
         r.squared =  map_dbl(stat, ~mcfadden(.)))%>%
  select(sumr, stat,r.squared)%>%
  unnest()%>%
  select(c(term:p.value, r.squared))

re_clutch%>%
  ggplot(aes(Female.SVL, rcm))+
  geom_point(aes(col = clutch.size), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed")+
  scale_y_continuous(name = "Relative clutch mass")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Clutch size")+
  theme(legend.text = element_text(size = 10))

re_clutch%>%
  select(Female.SVL, rcm)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(rcm ~ Female.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)

# calculating reproductive effort per embruo

re_embryo <- re%>%
  select(Field.Code, Embryo.Code, Egg.Weigth..g..,  Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
  group_by(Field.Code)%>%
  mutate(clutch.wt = sum(Egg.Weigth..g..), 
         Female.wt = Mother.Wt - clutch.wt)%>%
  group_by(Field.Code, Embryo.Code)%>%
  summarise(inv = Egg.Weigth..g../Female.wt,# investment per embryo
            Female.SVL = last(Female.SVL),
            Embryo.Sex = Embryo.Sex,
            Embryo.SVL = Embryo.SVL)%>%
  ungroup()

# using residuals to control for embryo development

emsvlinv <- lm(inv ~ Embryo.SVL, data = re_embryo)

# Plotting

re_embryo%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")

re_embryo%>%
  select(Female.SVL, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(inv ~ Female.SVL + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:p.value, r.squared)

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")+
  scale_shape_discrete(name = "Embryo Sex")

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  select(Embryo.Sex, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:p.value, r.squared)